#!/usr/bin/python
# -*- coding: utf-8 -*-

from numpy import matrix
from numpy.linalg import norm
import random
from itertools import izip

import psyco
psyco.full()

import os,ctypes
libKMeans=ctypes.CDLL(os.path.abspath(os.path.dirname(__file__))+"/KMeans.so")

class KMeansNode:
    def __init__(self,id,v,data):
        self.id=id
        self.v=matrix(v)
        self.data=data

class _cKMeansNode(ctypes.Structure):
    _fields_=[
        ("id",ctypes.c_int)
        ,("v",ctypes.POINTER(ctypes.c_float))
    ]

class _cKMeansCluster(ctypes.Structure):
    _fields_=[
        ("num",ctypes.c_int)
        ,("nodeIDs",ctypes.POINTER(ctypes.c_int))
    ]    
def cKMeansNode(id,v):
    cFloats=ctypes.c_float*len(v)
    vv=cFloats()
    for i in xrange(len(v)):
        vv[i]=float(v[i])
    node=_cKMeansNode(int(id),vv)
    node.dimension=len(v)
    return node

class cKMeansClustering:
    def __init__(self,nodeArray):
        self.nodes=nodeArray
    def Calc(self,teacherNodeArray):
        cKMeansClusteringFunc=libKMeans.KMeansClustering
        cKMeansClusteringFunc.restype=ctypes.POINTER(ctypes.POINTER(_cKMeansCluster))
        
        cfreeClusters=libKMeans.freeClusters;
        
        nodes=(ctypes.POINTER(_cKMeansNode)*len(self.nodes))()
        for i in xrange(len(self.nodes)):
            nodes[i]=ctypes.pointer(self.nodes[i])

        teachers=(ctypes.POINTER(_cKMeansNode)*len(teacherNodeArray))()
        for i in xrange(len(teacherNodeArray)):
            teachers[i]=ctypes.pointer(teacherNodeArray[i])
        retclusters=cKMeansClusteringFunc(nodes,ctypes.c_int(self.nodes[0].dimension),ctypes.c_int(len(self.nodes)),teachers,ctypes.c_int(len(teacherNodeArray)))
        
        ret=[[retclusters[k].contents.nodeIDs[i] for i in xrange(retclusters[k].contents.num)] for k in xrange(len(teacherNodeArray))]
        cfreeClusters(retclusters)
        return ret

    def Calc2(self,k):
        cKMeansClusteringFunc=libKMeans.KMeansPPClustering
        cKMeansClusteringFunc.restype=ctypes.POINTER(ctypes.POINTER(_cKMeansCluster))
        
        cfreeClusters=libKMeans.freeClusters;
        
        nodes=(ctypes.POINTER(_cKMeansNode)*len(self.nodes))()
        for i in xrange(len(self.nodes)):
            nodes[i]=ctypes.pointer(self.nodes[i])

        retclusters=cKMeansClusteringFunc(nodes,ctypes.c_int(self.nodes[0].dimension),ctypes.c_int(len(self.nodes)),ctypes.c_int(k))
        
        ret=[[retclusters[ki].contents.nodeIDs[i] for i in xrange(retclusters[ki].contents.num)] for ki in xrange(k)]
        cfreeClusters(retclusters)
        return ret
    
    def CalcX(self,k):
        cXMeansClusteringFunc=libKMeans.XMeansClustering
        cXMeansClusteringFunc.restype=ctypes.POINTER(ctypes.POINTER(_cKMeansCluster))
        
        cfreeClusters=libKMeans.freeClusters;
        
        nodes=(ctypes.POINTER(_cKMeansNode)*len(self.nodes))()
        for i in xrange(len(self.nodes)):
            nodes[i]=ctypes.pointer(self.nodes[i])

        pk=ctypes.pointer(ctypes.c_int(k))
        retclusters=cXMeansClusteringFunc(nodes,ctypes.c_int(self.nodes[0].dimension),ctypes.c_int(len(self.nodes)),pk)
        k=pk.contents.value
        
        ret=[[retclusters[ki].contents.nodeIDs[i] for i in xrange(retclusters[ki].contents.num)] for ki in xrange(k)]
        cfreeClusters(retclusters)
        return ret
class KMeansClustering:
    def __init__(self,nodeArray):
        self.nodes=nodeArray

    def __dist2Vec(self,i,j):
        if i==j: return 0
        if i>j: return self.__dist2Vec(j,i)
        if not self.__dCache.has_key(i): self.__dCache.setdefault(i,dict())
        if not self.__dCache[i].has_key(j): self.__dCache[i].setdefault(j,norm(self.__vecs[i]-self.__vecs[j]))
        
        return self.__dCache[i][j]
        
    def makeTeacherVectorArray(self,k):
        centers=[]
        self.__vecs=[nd.v for nd in self.nodes]
        self.__dCache=dict()
        
        # 1点目はランダム抽出
        idx=random.randint(0,len(self.__vecs)-1)
        centers.append(self.__vecs[idx]);

        #キョリ測
        dist=[self.__dist2Vec(i,idx) for i in xrange(len(self.__vecs))]

        for times in xrange(k-1):
            #$idx=argmin_{i}(sum_{j}(min(dist[j],|$vec[j]-$vec[i]|^2))
            vecsum=[sum((min(vecdists) for vecdists in izip(dist,(self.__dist2Vec(i,j) for j in xrange(len(self.__vecs)))))) for i in xrange(len(self.__vecs))]
            bestidx=vecsum.index(min(vecsum))

            centers.append(self.__vecs[bestidx])
            #キョリ測。もとより短ければ更新            
            #dist=[min(di) for di in izip(dist,(self.__dist2Vec(bestidx,i) for i in xrange(len(self.__vecs))))]
            for i in xrange(len(dist)):
                bestDist=self.__dist2Vec(bestidx,i)
                if dist[i]>bestDist: dist[i]=bestDist
            
        del self.__vecs
        del self.__dCache
        return centers
        
    def Calc2(self,k):
        '''K-Means++による初期値選択ののち、k個にクラスタリング'''
        
        #K-means++による初期値選択
        teachers=self.makeTeacherVectorArray(k)
        #クラスタリング
        return self.Calc(teachers)
    
    def Calc(self,teacherVectorArray):

        prevClusters=None
        while True:
            #教師の数だけクラスタを作る。
            clusters=[[] for t in teacherVectorArray]

            #最近傍のクラスタに入れる
            for node in self.nodes:
                nearestIdx=self.___calcNearestTeacherIndex(node,teacherVectorArray)
                clusters[nearestIdx].append(node)
            
            if not prevClusters is None:
                if self.__sameClusters(prevClusters,clusters):
                    #結果が同じなら終了。
                    del prevClusters
                    break
            
            #各クラスタの代表値を更新
            prevClusters=list(clusters)
            del clusters
            
            teacherVectorArray=[sum([node.v for node in cluster])/len(cluster) for cluster in prevClusters]
        
        #クラスタの代表値を格納
        self.teacherVectorResults=teacherVectorArray
        #結果を返す
        return clusters

    @staticmethod
    def ___calcNearestTeacherIndex(node,teacherVectorArray):
        """最近傍のクラスタのインデックスを求める"""        
        
        #クラスタの代表値との距離
        dist=[norm(node.v - teacher) for teacher in teacherVectorArray]
        return dist.index(min(dist))
    
    @staticmethod
    def __sameClusters(cs1,cs2):
        #for c1,c2 in izip(cs1,cs2):
            #print [nd.id for nd in c1]
            #print [nd.id for nd in c2]
        #print ""
        """クラスタ結果が同じかどうかを判別する"""
        for c1,c2 in izip(cs1,cs2):
            if len(c1)!=len(c2): return False
            
            #cluster1のノードそれぞれに対して
            for node1 in c1:
                #cluster2のどれかのノードとの対応が見つかればOK
                for node2 in c2:
                    if node1.id==node2.id: break
                else:
                    return False

        return True
    
    @staticmethod
    def _showclusters(clusters):
        for cluster in clusters:
            print [str(node.id) for node in cluster]


if __name__=="__main__":
    a=1

    nodes=[]
    if a==0:
        def p_u_s_h(v):
            nodes.append(KMeansNode(len(nodes),v,None))
    else:
        def p_u_s_h(v):
            nodes.append(cKMeansNode(1000+len(nodes),v))
        
    p_u_s_h((5.100, 3.500, 1.400, 0.200))
    p_u_s_h((4.900, 3.000, 1.400, 0.200))
    p_u_s_h((4.700, 3.200, 1.300, 0.200))
    p_u_s_h((4.600, 3.100, 1.500, 0.200))
    p_u_s_h((5.000, 3.600, 1.400, 0.200))
    p_u_s_h((5.400, 3.900, 1.700, 0.400))
    p_u_s_h((4.600, 3.400, 1.400, 0.300))
    p_u_s_h((5.000, 3.400, 1.500, 0.200))
    p_u_s_h((4.400, 2.900, 1.400, 0.200))
    p_u_s_h((4.900, 3.100, 1.500, 0.100))
    p_u_s_h((5.400, 3.700, 1.500, 0.200))
    p_u_s_h((4.800, 3.400, 1.600, 0.200))
    p_u_s_h((4.800, 3.000, 1.400, 0.100))
    p_u_s_h((4.300, 3.000, 1.100, 0.100))
    p_u_s_h((5.800, 4.000, 1.200, 0.200))
    p_u_s_h((5.700, 4.400, 1.500, 0.400))
    p_u_s_h((5.400, 3.900, 1.300, 0.400))
    p_u_s_h((5.100, 3.500, 1.400, 0.300))
    p_u_s_h((5.700, 3.800, 1.700, 0.300))
    p_u_s_h((5.100, 3.800, 1.500, 0.300))
    p_u_s_h((5.400, 3.400, 1.700, 0.200))
    p_u_s_h((5.100, 3.700, 1.500, 0.400))
    p_u_s_h((4.600, 3.600, 1.000, 0.200))
    p_u_s_h((5.100, 3.300, 1.700, 0.500))
    p_u_s_h((4.800, 3.400, 1.900, 0.200))
    p_u_s_h((5.000, 3.000, 1.600, 0.200))
    p_u_s_h((5.000, 3.400, 1.600, 0.400))
    p_u_s_h((5.200, 3.500, 1.500, 0.200))
    p_u_s_h((5.200, 3.400, 1.400, 0.200))
    p_u_s_h((4.700, 3.200, 1.600, 0.200))
    p_u_s_h((4.800, 3.100, 1.600, 0.200))
    p_u_s_h((5.400, 3.400, 1.500, 0.400))
    p_u_s_h((5.200, 4.100, 1.500, 0.100))
    p_u_s_h((5.500, 4.200, 1.400, 0.200))
    p_u_s_h((4.900, 3.100, 1.500, 0.200))
    p_u_s_h((5.000, 3.200, 1.200, 0.200))
    p_u_s_h((5.500, 3.500, 1.300, 0.200))
    p_u_s_h((4.900, 3.600, 1.400, 0.100))
    p_u_s_h((4.400, 3.000, 1.300, 0.200))
    p_u_s_h((5.100, 3.400, 1.500, 0.200))
    p_u_s_h((5.000, 3.500, 1.300, 0.300))
    p_u_s_h((4.500, 2.300, 1.300, 0.300))
    p_u_s_h((4.400, 3.200, 1.300, 0.200))
    p_u_s_h((5.000, 3.500, 1.600, 0.600))
    p_u_s_h((5.100, 3.800, 1.900, 0.400))
    p_u_s_h((4.800, 3.000, 1.400, 0.300))
    p_u_s_h((5.100, 3.800, 1.600, 0.200))
    p_u_s_h((4.600, 3.200, 1.400, 0.200))
    p_u_s_h((5.300, 3.700, 1.500, 0.200))
    p_u_s_h((5.000, 3.300, 1.400, 0.200))
    p_u_s_h((7.000, 3.200, 4.700, 1.400))
    p_u_s_h((6.400, 3.200, 4.500, 1.500))
    p_u_s_h((6.900, 3.100, 4.900, 1.500))
    p_u_s_h((5.500, 2.300, 4.000, 1.300))
    p_u_s_h((6.500, 2.800, 4.600, 1.500))
    p_u_s_h((5.700, 2.800, 4.500, 1.300))
    p_u_s_h((6.300, 3.300, 4.700, 1.600))
    p_u_s_h((4.900, 2.400, 3.300, 1.000))
    p_u_s_h((6.600, 2.900, 4.600, 1.300))
    p_u_s_h((5.200, 2.700, 3.900, 1.400))
    p_u_s_h((5.000, 2.000, 3.500, 1.000))
    p_u_s_h((5.900, 3.000, 4.200, 1.500))
    p_u_s_h((6.000, 2.200, 4.000, 1.000))
    p_u_s_h((6.100, 2.900, 4.700, 1.400))
    p_u_s_h((5.600, 2.900, 3.600, 1.300))
    p_u_s_h((6.700, 3.100, 4.400, 1.400))
    p_u_s_h((5.600, 3.000, 4.500, 1.500))
    p_u_s_h((5.800, 2.700, 4.100, 1.000))
    p_u_s_h((6.200, 2.200, 4.500, 1.500))
    p_u_s_h((5.600, 2.500, 3.900, 1.100))
    p_u_s_h((5.900, 3.200, 4.800, 1.800))
    p_u_s_h((6.100, 2.800, 4.000, 1.300))
    p_u_s_h((6.300, 2.500, 4.900, 1.500))
    p_u_s_h((6.100, 2.800, 4.700, 1.200))
    p_u_s_h((6.400, 2.900, 4.300, 1.300))
    p_u_s_h((6.600, 3.000, 4.400, 1.400))
    p_u_s_h((6.800, 2.800, 4.800, 1.400))
    p_u_s_h((6.700, 3.000, 5.000, 1.700))
    p_u_s_h((6.000, 2.900, 4.500, 1.500))
    p_u_s_h((5.700, 2.600, 3.500, 1.000))
    p_u_s_h((5.500, 2.400, 3.800, 1.100))
    p_u_s_h((5.500, 2.400, 3.700, 1.000))
    p_u_s_h((5.800, 2.700, 3.900, 1.200))
    p_u_s_h((6.000, 2.700, 5.100, 1.600))
    p_u_s_h((5.400, 3.000, 4.500, 1.500))
    p_u_s_h((6.000, 3.400, 4.500, 1.600))
    p_u_s_h((6.700, 3.100, 4.700, 1.500))
    p_u_s_h((6.300, 2.300, 4.400, 1.300))
    p_u_s_h((5.600, 3.000, 4.100, 1.300))
    p_u_s_h((5.500, 2.500, 4.000, 1.300))
    p_u_s_h((5.500, 2.600, 4.400, 1.200))
    p_u_s_h((6.100, 3.000, 4.600, 1.400))
    p_u_s_h((5.800, 2.600, 4.000, 1.200))
    p_u_s_h((5.000, 2.300, 3.300, 1.000))
    p_u_s_h((5.600, 2.700, 4.200, 1.300))
    p_u_s_h((5.700, 3.000, 4.200, 1.200))
    p_u_s_h((5.700, 2.900, 4.200, 1.300))
    p_u_s_h((6.200, 2.900, 4.300, 1.300))
    p_u_s_h((5.100, 2.500, 3.000, 1.100))
    p_u_s_h((5.700, 2.800, 4.100, 1.300))
    p_u_s_h((6.300, 3.300, 6.000, 2.500))
    p_u_s_h((5.800, 2.700, 5.100, 1.900))
    p_u_s_h((7.100, 3.000, 5.900, 2.100))
    p_u_s_h((6.300, 2.900, 5.600, 1.800))
    p_u_s_h((6.500, 3.000, 5.800, 2.200))
    p_u_s_h((7.600, 3.000, 6.600, 2.100))
    p_u_s_h((4.900, 2.500, 4.500, 1.700))
    p_u_s_h((7.300, 2.900, 6.300, 1.800))
    p_u_s_h((6.700, 2.500, 5.800, 1.800))
    p_u_s_h((7.200, 3.600, 6.100, 2.500))
    p_u_s_h((6.500, 3.200, 5.100, 2.000))
    p_u_s_h((6.400, 2.700, 5.300, 1.900))
    p_u_s_h((6.800, 3.000, 5.500, 2.100))
    p_u_s_h((5.700, 2.500, 5.000, 2.000))
    p_u_s_h((5.800, 2.800, 5.100, 2.400))
    p_u_s_h((6.400, 3.200, 5.300, 2.300))
    p_u_s_h((6.500, 3.000, 5.500, 1.800))
    p_u_s_h((7.700, 3.800, 6.700, 2.200))
    p_u_s_h((7.700, 2.600, 6.900, 2.300))
    p_u_s_h((6.000, 2.200, 5.000, 1.500))
    p_u_s_h((6.900, 3.200, 5.700, 2.300))
    p_u_s_h((5.600, 2.800, 4.900, 2.000))
    p_u_s_h((7.700, 2.800, 6.700, 2.000))
    p_u_s_h((6.300, 2.700, 4.900, 1.800))
    p_u_s_h((6.700, 3.300, 5.700, 2.100))
    p_u_s_h((7.200, 3.200, 6.000, 1.800))
    p_u_s_h((6.200, 2.800, 4.800, 1.800))
    p_u_s_h((6.100, 3.000, 4.900, 1.800))
    p_u_s_h((6.400, 2.800, 5.600, 2.100))
    p_u_s_h((7.200, 3.000, 5.800, 1.600))
    p_u_s_h((7.400, 2.800, 6.100, 1.900))
    p_u_s_h((7.900, 3.800, 6.400, 2.000))
    p_u_s_h((6.400, 2.800, 5.600, 2.200))
    p_u_s_h((6.300, 2.800, 5.100, 1.500))
    p_u_s_h((6.100, 2.600, 5.600, 1.400))
    p_u_s_h((7.700, 3.000, 6.100, 2.300))
    p_u_s_h((6.300, 3.400, 5.600, 2.400))
    p_u_s_h((6.400, 3.100, 5.500, 1.800))
    p_u_s_h((6.000, 3.000, 4.800, 1.800))
    p_u_s_h((6.900, 3.100, 5.400, 2.100))
    p_u_s_h((6.700, 3.100, 5.600, 2.400))
    p_u_s_h((6.900, 3.100, 5.100, 2.300))
    p_u_s_h((5.800, 2.700, 5.100, 1.900))
    p_u_s_h((6.800, 3.200, 5.900, 2.300))
    p_u_s_h((6.700, 3.300, 5.700, 2.500))
    p_u_s_h((6.700, 3.000, 5.200, 2.300))
    p_u_s_h((6.300, 2.500, 5.000, 1.900))
    p_u_s_h((6.500, 3.000, 5.200, 2.000))
    p_u_s_h((6.200, 3.400, 5.400, 2.300))
    p_u_s_h((5.900, 3.000, 5.100, 1.800))

    if a==0:
        #KMeansClustering._showclusters(KMeansClustering(nodes).Calc2(5))
        for r in KMeansClustering(nodes).Calc2(5):
            print [nd.id for nd in r]
    else:
        ret=cKMeansClustering(nodes).CalcX(10)
        ret.sort(key=lambda r:r[0])
        for r in ret:
            r.sort()
            print r